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Abstract 

The energy landscapes of electrostatically charged particles embedded on constant mean curvature sur- 
faces are analysed for a wide range of system size, curvature, and interaction potentials. The surfaces are 
taken to be rigid, and the basin-hopping method is used to locate the putative global minimum structures. 
The defect motifs favoured by potential energy agree with experimental observations for colloidal systems: 
extended defects (scars and pleats) for weakly positive and negative Gaussian curvatures, and isolated de- 
fects for strongly negative Gaussian curvatures. Near the phase boundary between these regimes the two 
motifs are in strong competition, as evidenced from the appearance of distinct funnels in the potential energy 
landscape. We also report a novel defect motif consisting of pentagon pairs. 



1 



Introduction - The distribution of electrostatically charged particles on curved surfaces [1] pro- 
vides a testing ground for global optimisation algorithms D2k7|], as well as useful insights into a 



number of materials science and biological applications, inc 



%, fullerene structures 



JlO, 



uding the packing of virus capsids 



llll . colloidal crystals [|12l 1 1 3fl . and proteins on lipid membranes 



1411 . A key issue in understanding the resulting structures is the appearance of defects, which 



strongly influence the mechanical, electronic and optical properties. Defect structures have even 



been exploited for rational design of templated se 



been used to create divalent metal nanoparticles [15] 



f-assembly; for example, polar defects have 



Frustration is intrinsic to curved surfaces because of the competition between local order and 
long-range geometrical constraints, and defects appear to screen the resulting geometrical stresses, 
even for the ground state configuration. For example, while a hexagonal lattice is possible for 
a flat surface, this arrangement cannot occur on a sphere without introducing isolated five-fold 
disclinations (pentagons; with positive topological charges), or composite structures with seven- 

170 . Generally, adjacent 



fold disclinations (heptagons; negative topological charges) Q2H7L I 
combinations of pentagons and heptagons may appear as topologically uncharged or charged lines 
of dislocations, corresponding to pleats (e.g. pentagon-heptagon topological dipoles) and scars 
(e.g. pentagon-heptagon-pentagon arrangements), respectively. 

Isolated heptagons may also exist if the Gaussian curvature of the substrate is strongly negative, 
as was shown recently in studies of two-dimensional colloidal crystals on the surface of capillary 
bridges Jl8]. By systematically varying the shape and thus the curvature of the substrate, a se- 
quence of transitions was observed from zero defects to isolated dislocations, pleats, scars, and 
isolated disclinations. 

In the present contribution we characterise the energy landscape of electrostatically charged 
particles embedded on constant mean curvature surfaces, modelling in particular the recent ex- 
perimental setup by Irvine et al. for colloidal crystals 11181. LL9Q . Employing the basin-hopping 
algorithm [|20l- l23|l . we identify likely global minimum configurations for a wide range of system 
sizes and surface curvatures, considering several different interaction potentials. Not only are we 
able to reproduce the experimental sequence of defect transitions, but we also identify a new de- 
fect motif corresponding to pentagon pairs, which may appear on surfaces with both positive and 
negative Gaussian curvatures. Furthermore, we show here for the first time the hierarchical na- 
ture of the potential energy landscape for these systems. Using disconnectivity graphs Q24[], we 
demonstrate the appearance of separate funnels in the landscape, corresponding to competing de- 
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feet motifs favouring dislocations and disclinations. Our analysis also provides information on the 
rearrangement mechanisms between defect patterns, as well as insight into likely thermodynamic 
signatures for structural transitions. 

Methodology - We consider N identical electrostatically charged particles embedded on the 
surface of catenoids (zero mean curvature) and unduloids (non-zero mean curvature) [|25ll . To rep- 

n 

resent the interactions of these colloids trapped on fluid interfaces H26H . we have mainly considered 
Yukawa potentials of the form: V = J2f>i T~ e~( rij / x \ where is the Euclidean distance be- 
tween particle i and j, and A is the screening length. We choose A = 0.1 so that the ratio between 
the screening length and the length of the capillary bridges is of the same order as that used in 



experiments [|18 



26D. 



We have also considered four other potentials: Yukawa with A = 0.4, Coulomb, Lennard- 
Jones, and a repulsive Lennard- Jones form. Qualitatively, the defect motifs are the same provided 
that the particle density per unit area is sufficiently high, indicating that our defect analysis should 
have wide-ranging applicability to curved surfaces. This regime corresponds to N > 200. If the 
number of electrostatically charged particles is smaller, they tend to populate the boundaries when 
Coulomb or other similarly long-range potentials are used; for the Lennard- Jones potential, we 
observe patches of hexagonal lattice separated by large spaces. We therefore focus on 200 < N < 
600 in the present work. 

Any point on the surface of an unduloid or a catenoid can be parameterised using two 
variables u and v I127n . For catenoids, the mapping to Cartesian coordinates takes the fol- 
lowing form (x, y,z) = (ccosh (t>/c) cos (u), c cosh (v/c) svs\{u),v), where < u < 2ix, 
—z m < v < z m , and c is a free parameter corresponding to the waist of the catenoid in the 
z = plane. The corresponding transformation is more elaborate for unduloids (x, y, z) = 
((m sin fiv + n) 1 / 2 cos (u), [m sin [iv + n) 1 / 2 sin (u), aF (fJ,v/2 — 7r/4, k) + cE (fJ,v/2 — 7r/4, k)), 
with n = 2/(a + c), k 2 = (c 2 - a 2 )/c 2 , m = (c 2 - a 2 )/2, and n = (c 2 + a 2 )/2. k) and 
E((j), k) are elliptic integrals of the first and second kind. We tune the shape of the unduloids by 
varying the parameters a, c, and the range of values for v Q2VU - For a neck shape unduloid [e.g. 
Fig. dta-c)], v is centred around v c = 37r/2/x, while for a barrel shape unduloid, v c = 7r/2/i. The 
maximum and minimum values of v are chosen so that — z m < z < z rn . 

To characterise the most favourable geometries we employ basin-hopping global optimisation 
II20h23I1 . In this method, random geometrical perturbations are followed by energy minimisation, 
and moves are accepted or rejected based upon the energy differences between local minima. The 



3 



minimisation procedure transforms the energy landscape of the system into the set of catchment 
basins for the local minima, and makes the basin-hopping method very effective for finding low- 
lying structures. Further details are provided in the supporting information. 

We have also constructed databases of connected minima for selected systems, starting from the 
low-lying structures found in the basin-hopping runs. We employed doubly-nudged elastic band 
transition state searches [|28ll . where images corresponding to local maxima were tightly converged 
to transition states using hybrid eigenvector-following [2^, 3(3]. This procedure provides a global 
survey of the potential energy landscape, which we then visualise using disconnectivity graphs 
112411 . We find that distinct structural arrangements of the particles result in separate funnels in the 
landscape. Moreover, the database of connected minima and transition states allows us to predict 
energy barriers and rearrangement mechanisms between defect rearrangements. 

To visualise the defect structures, we use Voronoi constructions. Thus, pentagons, hexagons, 
and heptagons correspond to particles with five, six, and seven neighbours. All the results pre- 
sented here were obtained using the GMIN, OPTIM, and PATHS AMPLE programs Q3JJ] , which 
are available for use under the GNU public license. 

Defect motifs on unduloids - We first analyse the arrangement of electrostatically charged parti- 
cles embedded on non-zero constant mean curvature surfaces (unduloids). The results for N = 600 
identical particles are presented in Fig. [TJ The unduloid parameters a and c were varied while keep- 
ing the height (2 z m = 1.5) and volume (V = 2.65) of the unduloids constant. Experimentally, 



220. 



these parameters correspond to capillary bridges for surfaces with different contact angles [|32L 
The defect motifs and sequence of transitions were found to be consistent as we varied the number 
of particles, and specific results for N = 200 are presented in the supporting information. 

Due to the finite number of particles edge effects are an intrinsic feature, and we find that 
the first two rows of Voronoi cells from the boundary consistently have more defects. As for 
the experimental results on colloidal crystals fl 1 811 we find a series of defect transitions, including 
isolated heptagons (Fig.[TJ). The Gaussian curvature is not constant on the surface of an unduloid, 
and it is most negative at the waist, where isolated heptagons are located. Far from the waist and 
near the edge of the unduloid, we observe lines of dislocations. 

For less negative total Gaussian curvatures [Figs Qlb-c)], isolated heptagons disappear and 
dislocation lines prevail. The length of the dislocations is also found to correlate strongly with the 
curvature: the length is shorter when the curvature is less negative. For weakly negative Gaussian 
curvatures, the pentagon-heptagon dipole is a common motif, together with a pair of pentagons. 
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(a) [0.50, 885.89] 



(b) [0.60, 8.75] 



(c) [0.70, 1.57] 



(d) [0.75, 0.75] 



(e) [0.06, 0.85] 



Negative Gaussian Curvature 
(f) ^'♦'"'''ilTST 



Positive Gaussian Curvature 




FIG. 1: Defect motifs for non-zero constant mean curvature surfaces, (a-e) The number of electrostatically 
charged particles is N = 600, and the curvature of the capillary bridges is varied while keeping the height 
and volume of the bridges constant. The corresponding values for the unduloid parameters [a, c] are given in 
square brackets, (f-i) Comparisons to experimental results of Irvine et al. for colloidal crystals, reproduced 
(will seek permission) from 111811 . Here, red particles have seven neighbors, while yellow particles have five 
neighbors. In the experiments, the curvature is tuned by stretching/compressing the liquid bridge, and the 
sequence of defect transitions observed is identical to our global optimisation results. 

The pentagon pair motif is surrounded by seven hexagons, as for an isolated heptagon, and to the 
best of our knowledge, this motif has not been reported before. Since isolated pentagons have a 
positive topological charge, it is somewhat surprising that double pentagons may exist in a bound 
state, and that they can be favourable for surfaces with negative Gaussian curvature. 

We observe no defects for cylinders [zero Gaussian curvature, Fig.QJd)], as expected. Dislo- 
cation lines as well as pentagon pairs then reappear for unduloid with a positive total Gaussian 
curvature [Fig. (He)]. The orientations of the dislocations and pentagon pairs are reversed for sur- 
faces with positive and negative Gaussian curvature. As predicted by continuum elastic theories 



34f | . for negative curvatures, the 7-5 dislocation axis runs along the meridian, pointing in 



the direction of the boundary. For positive Gaussian curvatures, the dislocation axis points to the 
centre of the unduloid instead. Similarly, the pentagon pairs point to the boundary for negative 
Gaussian curvatures, and to the centre of the unduloids for positive Gaussian curvatures. 

Defect motifs on catenoids - We have also analysed the defect structures for a family of 
catenoids (minimal surfaces) with varying values of the waist parameter c. We adjust the height 
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(a) c = 0.4 




(b) c = 0.5 



m 



(d) c = 0.7 



(e)c = 0.8 




Increasingly negative total Gaussian curvature 



FIG. 2: Defect motifs on catenoid surfaces for N = 600, with variable waist radius c. 



parameter z m so that the radius of the catenoid at z m , ccosh (z m /c) = 1. This procedure results in 
less negative total Gaussian curvatures with increasing waist radius c, see Fig. |2] 

For strongly negative Gaussian curvature (small c in Fig. |2]), we always observe isolated hep- 
tagons at the waist of the catenoids, independent of N. We illustrate the disconnectivity graph 
113511 for N = 600 and c = 0.40 in Fig. Ufa). The potential energy landscape is clearly hierarchical 
in character. Interestingly, there are only two dominant defect configurations near the waist with 
a high degree of symmetry, which we label as a and (3 in the inset, each consisting of eight hep- 
tagons. The main variation in the structural arrangements of the ions comes from the dislocation 
lines near the edges (see the Supporting Information for animations). The energy barrier for in- 
terconverting waist configurations a and (3 is approximately AEi ~ 0.4 (in reduced units) if the 
configuration near the edge of the catenoid is roughly preserved. On the other hand, the barrier 
for rearranging the particle distribution near the edges can be much higher, AE 2 ~ 2.2 (reduced 
units), as indicated in Fig. Ufa). From the pathway we see that the high barrier is due to global 
concerted rotation of several layers of ions near the edges of the catenoids. 

For weakly negative Gaussian curvatures (large c in Fig. [2]), we never find isolated heptagons, 
and only lines of dislocations and pentagon pairs exist. It is worth noting that these defect motifs 
have the same orientations as for unduloids with negative Gaussian curvature, above. The energy 
landscape for this parameter regime is also simpler. In particular, we find that the number of 
possible local minima is considerably smaller. 

The transition from defect patterns favouring disclinations to dislocations occurs at moderately 
negative Gaussian curvatures, c ~ 0.55, which we determine by constructing a defect phase dia- 
gram as a function of the number of electrostatically charged particles, N, and the catenoid waist 
radius, c (Fig. 0). Our global optimisation results further suggest that this transition is only weakly 
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FIG. 3: Defect phase diagram as a function of the catenoid waist radius, c, and number of electrostatically 
charged particles, N. Squares denote global minimum structures containing isolated heptagons, while 
circles correspond to structures without isolated heptagons. 

dependant on the number of particles in this size range. The phase diagram for the repulsive 
Lennard- Jones potential is shown in the Supporting Information, where we find similar behaviour, 
except that the phase boundary is shifted to larger c. Our results for N < 600 are probably far 
from the continuum limit, where Bowick and Yao Il34ll predict that c will decrease monotonically 
with N. Additionally, the phase boundary found in Fig. [3] is neither smooth nor monotonic. We 
can explain this observation by analysing the disconnectivity graph for N = 600 and c = 0.575. 
As shown in Fig. BJb), there are two distinct funnels, one favouring dislocations and the other 
favouring disclinations. The energy difference between the lowest configurations in the two fun- 
nels is very small compared to the energy barrier (AE 3 ~ 0.7 in reduced units) associated with 
interconverting the two defect motifs. 

Discussion - We have investigated the most favourable defect motifs for electrostatically 
charged particles embedded on zero and non-zero constant mean curvature surfaces. By vary- 
ing the curvature of the embedding surface, we have characterised the detailed structures and 
energetics of a wide range of defect motifs, including dislocation lines and isolated disclinations, 
in excellent agreement with experiment and predictions from continuum elastic theories. The 
appearance of a new defect motif consisting of pentagon pairs is also predicted. 

Taking the typical experimental values for electrostatic interactions between colloids embedded 
on interfaces, the energy scale for one simulation unit corresponds to between 10 and 100k B T (at 
room temperature). The energy barriers of interest, as shown in Fig. H are therefore quite large 
and the structures observed in experiments could be trapped in local minima. For experimen- 
tally relevant temperatures the favoured defect motifs are primarily determined by the potential 
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(a) 



(b) 



M 




FIG. 4: Disconnectivity graphs for N = 600 electrostatically charged particles embedded on catenoids with 
waist radius (a) c = 0.40 and (b) c = 0.575. The insets in (a) and (b) show representative minima for the 
corresponding funnels in the potential energy landscape. We also indicate the typical energy barriers for 
interconverting different defect motifs. AEi, AE2, and AE3 are respectively 0.4, 2.2, and 0.7 in reduced 
units. 



energy, which is consistent with the successful structure predictions that we have reported. We 
hope these results will stimulate future theoretical and experimental work. In particular, it will be 
interesting to modulate (reduce) the strength of interactions between the electrostatically charged 
colloids, so that the defect rearrangement mechanisms and thermodynamics of the system may 
be fully explored. Since the energy landscape is hierarchical, with different energy scales for 
interconverting the defect motifs, we expect the thermodynamics of the system to be very rich. 
For example, we predict that there will be multiple signatures in the heat capacity corresponding 
to alternative defect morphologies D23|]. Calculations of the thermodynamics, kinetic rates, and 
rearrangement pathways for these defect motif transitions are currently underway and will be pre- 
sented elsewhere. The interplay between the arrangement of the electrostatically charged particles 
and possible deformation of the interface is still an open question, and allowing for flexible rather 
than rigid curved surfaces is another avenue for future research. 



Analog models of the defect motifs presented here can be constructed using Polydron tiles [36], 



as illustrated in the Supporting Information (SI). The SI also contains animations of the defect 
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morphologies shown in Figs[T]|2]and|H me tabulated numbers of positive and negative topological 
charges, and detailed analysis of the net topological charges as a function of the integrated Gaus- 
sian curvature for several representative cases. These putative global minima found here will be 
made available online from the Cambridge Cluster Database Q3VU - 
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